Diferenciación numérica

En esta práctica exploraremos las formulas para la diferenciación numérica, así como un par de maneras de hacer los calculos de estas mas eficientes.

Empecemos con una función compleja de derivar:

$$ f(x) = \tanh{\left(\ln{\left(x^x\right)}\right)} $$

estas funciones se encuentran en la liberria de numpy:


In [ ]:
from numpy import tanh, log, linspace

y definimos la función de manera simple con la notacion lambda:


In [ ]:
f = lambda x: tanh(log(x**x))

De manera preliminar, podemos graficar unos cuantos puntos para darnos una idea de la forma general de la función, para esto definimos unos cuantos puntos:


In [ ]:
datos_x = linspace(0.2, 1.8, 9)
datos_x

y obtenemos el valor de estos datos en la función a derivar:


In [ ]:
datos_y = f(datos_x)

In [ ]:
from matplotlib.pyplot import plot
%matplotlib inline

In [ ]:
plot(datos_x, datos_y)

Por la gráfica podemos ver que para el valor de $x=1$, la grafica pasa aproximadamente por $0$, evaluando en este punto, vemos que realmente es:


In [ ]:
f(1)

Pero mas importante aun, este es un ejercicio del libro de texto, por lo que sabemos que en ese punto la derivada de esta función tiene un valor $f'(1) = 1$, por lo que estaremos viendo aproximaciones de este valor.


In [ ]:
plot(datos_x, datos_y)
plot([0.2,1.8], [-0.8,0.8])

Una vez que tenemos claro nuestro objetivo, empezamos a escribir las funciones que calcularán el valor de la derivada en cierto punto, por ejemplo, podemos empezar con las funciones que calculan las derivadas de dos puntos:

$$ f'(x_0) = \frac{f(x_1) - f(x_ 0)}{h} $$$$ f'(x_0) = \frac{f(x_0) - f(x_{-1})}{h} $$

en donde el valor de $h$ es calculado con el paso de los datos, o bien la diferencia entre $x_0$ y $x_1$, por facilidad de implementación, estas las calcularemos con $h= \left| x_1 - x_0 \right|$

Note que el conjunto de valores que hemos utilizado, tiene un valor constante de $h$, es importante que los datos en diferenciación numérica sean siempre calculados asi.


In [ ]:
def derivada_adelante_dos(func, x0, x1):
    return (func(x1) - func(x0))/abs(x1-x0)

In [ ]:
derivada_adelante_dos(f, 1, 1.2)

In [ ]:
def derivada_atras_dos(func, x0, x1):
    return (func(x0) - func(x1))/abs(x1-x0)

In [ ]:
derivada_atras_dos(f, 1, 0.8)

Lo mas evidente de estos resultados es que no son exactos, mas aún, ni siquiera estan perfectamente centrados en el valor principal, podemos ver esto aun mas evidente con las formulas para la derivación con 3 puntos:

$$ f'(x_0) = \frac{-3f(x_0) + 4f(x_1) - f(x_2)}{2h} $$$$ f'(x_0) = \frac{3f(x_0) - 4f(x_{-1}) + f(x_{-2})}{2h} $$

In [ ]:
def derivada_adelante_tres(func, x0, x1, x2):
    return (-3*func(x0) + 4*func(x1) - f(x2))/abs(x2-x0)

In [ ]:
derivada_adelante_tres(f, 1, 1.2, 1.4)

In [ ]:
def derivada_atras_tres(func, x0, x1, x2):
    return (3*func(x0) - 4*func(x1) + f(x2))/abs(x2-x0)

In [ ]:
derivada_atras_tres(f, 1, 0.8, 0.6)

Por lo que si quisieramos calcular el error, restando el valor de uno con el otro, nos hubiera reportado un error menor al realmente ocurrido.

Con esto podemos concliur que no hay una manera analítica de determinar el error obtenido, tan solo tenemos la garantía de que las formulas de 2, 3 y 5 puntos nos van a entregar errores del orden de magnitud $\mathcal{O}\left( h \right)$, $\mathcal{O}\left( h^2 \right)$ y $\mathcal{O}\left( h^4 \right)$; es decir, al menos tenemos la garantia de que usar las formulas de 5 puntos nos entregarán un menor error, siempre y cuando $h<1$.

Definamos pues las formulas para las derivadas de 5 puntos:

$$ f'(x_0) = \frac{-25f(x_0) + 48f(x_1) - 36f(x_2) + 16f(x_3) - 3f(x_4)}{12h} $$$$ f'(x_0) = \frac{25f(x_0) - 48f(x_{-1}) + 36f(x_{-2}) - 16f(x_{-3}) + 3f(x_{-4})}{12h} $$

In [ ]:
def derivada_adelante_cinco(func, x0, x1, x2, x3, x4):
    return (-25*func(x0) + 48*func(x1) - 36*func(x2) + 16*func(x3) - 3*func(x4))/(3*abs(x4 - x0))

In [ ]:
derivada_adelante_cinco(f, 1, 1.2, 1.4, 1.6, 1.8)

In [ ]:
def derivada_atras_cinco(func, x0, x1, x2, x3, x4):
    return (25*func(x0) - 48*func(x1) + 36*func(x2) - 16*func(x3) + 3*func(x4))/(3*abs(x4 - x0))

In [ ]:
derivada_atras_cinco(f, 1, 0.8, 0.6, 0.4, 0.2)

Sin embargo, estas formulas utilizan una y otra vez la funcion original... y esto es adecuado si no te preocupa el poder de procesamiento del dispositivo en el que se calcula, esto es si tienes acceso a la función analítica.

En este ejemplo teniamos una formula para la función a derivar, sin embargo va a haber ocasiones en las que queremos saber la derivada de un conjunto de datos, en este caso, no sabremos la forma analítica de la función original y por lo tanto, este enfoque es inutil.

Introduciremos el concepto de memoización, el cual consiste en calcular todo lo posible de antemano y utilizarlo cuando nos sea necesario.

En este caso, ya tenemos un conjunto de datos que calculamos anteriormente, datos_x y datos_y fueron calculados tomando en cuenta $h$ y $f(x)$, por lo que realmente no necesitamos estos valores para nuestras funciones.

Si definimos las funciones para las derivadas con dos puntos, tendremos lo siguiente:


In [ ]:
def dadel2(xs, ys, i):
    return (ys[i+1] - ys[i])/abs(xs[i+1]-xs[i])

In [ ]:
dadel2(datos_x, datos_y, 4)

In [ ]:
def datr2(xs, ys, i):
    return (ys[i] - ys[i-1])/abs(xs[i] - xs[i-1])

In [ ]:
datr2(datos_x, datos_y, 4)

y por los resultados que nos dan, podemos ver que son equivalentes.

Problemas

  1. Obtenga funciones para la obtencion de las derivadas numéricas de 3 y 5 puntos, hacia adelante, centrales y hacia atras, utilizando la técnica de memoización (utiliza las formulas ubicadas en las páginas 456 y 457 de tu libro de texto.
  2. Resuelva el ejercicio 5.41 del libro de texto.
  3. Utilizando el siguiente código para obtener el conjunto de datos, obtenga la derivada de la señal correspondiente a la posición de un servomotor para obtener su velocidad.

In [ ]:
import json
datos = None
with open('datos.json') as archivo:
     datos = json.load(archivo)
        
posicion, tiempo = datos[0], datos[1]